! Double precision complex error function
! by Alan.Miller @ vic.cmis.csiro.au
! http://www.ozemail.com.au/~milleraj

pure function zerf(z,mod) result(w)
  complex(8),intent(in) :: z
  complex(8)            :: w
  real(8),dimension(2)  :: z_,w_
  integer,optional,intent(in)      :: mod
  integer               :: mod_
  mod_=0;if(present(mod))mod_=mod
  z_(1)=real(z,8)
  z_(2)=aimag(z)
  call dcerf(mod_,z_,w_)
  w=cmplx(w_(1),w_(2))
end function zerf



PURE SUBROUTINE dcerf (mo, z, w)
  !-----------------------------------------------------------------------
  !               COMPUTATION OF THE COMPLEX ERROR FUNCTION
  !                          -----------------
  !                       W = ERF(Z)    IF MO = 0
  !                       W = ERFC(Z)   OTHERWISE
  !-----------------------------------------------------------------------
  INTEGER, INTENT(IN)     :: mo
  REAL (8), INTENT(IN)   :: z(2)
  REAL (8), INTENT(OUT)  :: w(2)
  REAL (8) :: m, n, n2, n4, np1
  REAL (8) :: c2, d, d2, e, eps, r, sn, tol, x, y
  REAL (8) :: a0(2), an(2), b0(2), bn(2)
  REAL (8) :: g0(2), gn(2), h0(2), hn(2)
  REAL (8) :: qf(2), sm(2), sz(2), tm(2), ts(2), w0(2), wn(2)
  !------------------------
  !     C = 1/SQRT(PI)
  !------------------------
  REAL (8), PARAMETER :: c = .56418958354775628694807945156077d0
  !------------------------
  !     ****** EPS IS A MACHINE DEPENDENT CONSTANT. EPS IS THE
  !            SMALLEST NUMBER SUCH THAT 1.d0 + EPS > 1.d0.
  eps = EPSILON(1.0d0)
  !------------------------
  x = z(1)
  y = z(2)
  sn = 1.d0
  IF (x >= 0.d0) GO TO 10
  x = -x
  y = -y
  sn = -1.d0
10 r = x*x + y*y
  sz(1) = x*x - y*y
  sz(2) = 2.d0*x*y
  IF (r <= 1.d0) GO TO 20
  IF (r >= 144.d0) GO TO 100
  IF (ABS(y) > 31.8d0*x) GO TO 50
  IF (ABS(y) > 7.0d0*x .AND. r < 64.d0) GO TO 50
  IF (ABS(y) > 3.2d0*x .AND. r < 49.d0) GO TO 50
  IF (ABS(y) > 2.0d0*x .AND. r < 36.d0) GO TO 50
  IF (ABS(y) > 1.2d0*x .AND. r < 25.d0) GO TO 50
  IF (ABS(y) > 0.9d0*x .AND. r < 16.d0) GO TO 50
  IF (r >= 6.25d0) GO TO 80
  IF (ABS(y) > 0.6d0*x) GO TO 50
  IF (r >= 4.0d0) GO TO 40
  d = x - 2.d0
  IF (d*d + y*y < 1.d0) GO TO 40
  GO TO 50
  !                          TAYLOR SERIES
20 c2 = c + c
  tm(1) = c2*x
  tm(2) = c2*y
  sm(1) = tm(1)
  sm(2) = tm(2)
  tol = 2.d0*eps
  m = 0.d0
21 m = m + 1.d0
  d = m + m + 1.d0
  ts(1) = tm(1)*sz(1) - tm(2)*sz(2)
  ts(2) = tm(1)*sz(2) + tm(2)*sz(1)
  tm(1) = -ts(1)/m
  tm(2) = -ts(2)/m
  ts(1) = tm(1)/d
  ts(2) = tm(2)/d
  sm(1) = sm(1) + ts(1)
  sm(2) = sm(2) + ts(2)
  IF (dnorm(ts(1),ts(2)) > tol*dnorm(sm(1),sm(2))) GO TO 21
  IF (mo /= 0) GO TO 30
  w(1) = sn*sm(1)
  w(2) = sn*sm(2)
  RETURN
30 IF (sn == 1.d0) GO TO 31
  w(1) = 1.d0 + sm(1)
  w(2) = sm(2)
  RETURN
31 w(1) = 0.5d0 + (0.5d0 - sm(1))
  w(2) = -sm(2)
  RETURN
  !                  TAYLOR SERIES AROUND Z0 = 2
40 tm(1) = x
  tm(2) = y
  CALL erfcm2 (0, tm, w)
  IF (mo /= 0) GO TO 41
  w(1) = sn*(0.5d0 + (0.5d0 - w(1)))
  w(2) = - sn*w(2)
  RETURN
41 IF (sn > 0.d0) RETURN
  w(1) = 2.d0 - w(1)
  w(2) = - w(2)
  RETURN
  !            PADE APPROXIMATION FOR THE TAYLOR SERIES
  !                    FOR  (EXP(Z*Z)/Z)*ERF(Z)
50 d = 4.d0
  IF (r > 16.d0) d = 16.d0
  IF (r > 64.d0) d = 64.d0
  d2 = d*d
  CALL dcrec (sz(1), sz(2), w(1), w(2))
  a0(1) = 1.d0
  a0(2) = 0.d0
  an(1) = (w(1) + 4.d0/15.d0)*d
  an(2) = w(2)*d
  b0(1) = 1.d0
  b0(2) = 0.d0
  bn(1) = (w(1) - 0.4d0)*d
  bn(2) = w(2)*d
  CALL cdivid (an(1), an(2), bn(1), bn(2), wn(1), wn(2))
  tol = 10.d0*eps
  n4 = 0.d0
60 n4 = n4 + 4.d0
  e = (n4 + 1.d0)*(n4 + 5.d0)
  tm(1) = d*(w(1) - 2.d0/e)
  tm(2) = d*w(2)
  e = d2*(n4*(n4 + 2.0))/((n4 - 1.0)*(n4 + 3.0)*(n4 + 1.0)**2)
  qf(1) = (tm(1)*an(1) - tm(2)*an(2)) + e*a0(1)
  qf(2) = (tm(1)*an(2) + tm(2)*an(1)) + e*a0(2)
  a0(1) = an(1)
  a0(2) = an(2)
  an(1) = qf(1)
  an(2) = qf(2)
  qf(1) = (tm(1)*bn(1) - tm(2)*bn(2)) + e*b0(1)
  qf(2) = (tm(1)*bn(2) + tm(2)*bn(1)) + e*b0(2)
  b0(1) = bn(1)
  b0(2) = bn(2)
  bn(1) = qf(1)
  bn(2) = qf(2)
  w0(1) = wn(1)
  w0(2) = wn(2)
  CALL cdivid (an(1), an(2), bn(1), bn(2), wn(1), wn(2))
  IF (dnorm(wn(1) - w0(1), wn(2) - w0(2)) > tol*dnorm(wn(1), wn(2))) GO TO 60
  c2 = c + c
  sm(1) = c2*(x*wn(1) - y*wn(2))
  sm(2) = c2*(x*wn(2) + y*wn(1))
  e = EXP(-sz(1))
  qf(1) = e*COS(-sz(2))
  qf(2) = e*SIN(-sz(2))
  tm(1) = qf(1)*sm(1) - qf(2)*sm(2)
  tm(2) = qf(1)*sm(2) + qf(2)*sm(1)
  w(1) = sn*tm(1)
  w(2) = sn*tm(2)
  IF (mo == 0) RETURN
  w(1) = 1.d0 - w(1)
  w(2) = - w(2)
  RETURN
  !         PADE APPROXIMATION FOR THE ASYMPTOTIC EXPANSION
  !                    FOR  Z*EXP(Z*Z)*ERFC(Z)
80 d = 4.d0*r
  IF (r < 16.d0) d = 16.d0*r
  d2 = d*d
  tm(1) = sz(1) + sz(1)
  tm(2) = sz(2) + sz(2)
  g0(1) = 1.d0
  g0(2) = 0.d0
  gn(1) = (2.d0 + tm(1))/d
  gn(2) = tm(2)/d
  h0(1) = 1.d0
  h0(2) = 0.d0
  tm(1) = 3.d0 + tm(1)
  hn(1) = tm(1)/d
  hn(2) = tm(2)/d
  CALL cdivid (gn(1), gn(2), hn(1), hn(2), wn(1), wn(2))
  np1 = 1.d0
  tol = 10.d0*eps
90 n = np1
  np1 = n + 1.d0
  n2 = n + n
  e = (n2*(n2 + 1.d0))/d2
  tm(1) = tm(1) + 4.d0
  qf(1) = (tm(1)*gn(1) - tm(2)*gn(2))/d - e*g0(1)
  qf(2) = (tm(1)*gn(2) + tm(2)*gn(1))/d - e*g0(2)
  g0(1) = gn(1)
  g0(2) = gn(2)
  gn(1) = qf(1)
  gn(2) = qf(2)
  qf(1) = (tm(1)*hn(1) - tm(2)*hn(2))/d - e*h0(1)
  qf(2) = (tm(1)*hn(2) + tm(2)*hn(1))/d - e*h0(2)
  h0(1) = hn(1)
  h0(2) = hn(2)
  hn(1) = qf(1)
  hn(2) = qf(2)
  w0(1) = wn(1)
  w0(2) = wn(2)
  CALL cdivid (gn(1), gn(2), hn(1), hn(2), wn(1), wn(2))
  IF (dnorm(wn(1) - w0(1), wn(2) - w0(2)) > tol*dnorm(wn(1), wn(2))) GO TO 90
  tm(1) = x*hn(1) - y*hn(2)
  tm(2) = x*hn(2) + y*hn(1)
  CALL cdivid (c*gn(1), c*gn(2), tm(1), tm(2), sm(1), sm(2))
  GO TO 130
  !                      ASYMPTOTIC EXPANSION
100 CALL dcrec (x, y, tm(1), tm(2))
  sm(1) = tm(1)
  sm(2) = tm(2)
  qf(1) = tm(1)*tm(1) - tm(2)*tm(2)
  qf(2) = 2.d0*tm(1)*tm(2)
  tol = 2.d0*eps
  d = -0.5d0
110 d = d + 1.d0
  ts(1) = tm(1)*qf(1) - tm(2)*qf(2)
  ts(2) = tm(1)*qf(2) + tm(2)*qf(1)
  tm(1) = -d*ts(1)
  tm(2) = -d*ts(2)
  sm(1) = sm(1) + tm(1)
  sm(2) = sm(2) + tm(2)
  IF (dnorm(tm(1), tm(2)) > tol*dnorm(sm(1), sm(2))) GO TO 110
  sm(1) = c*sm(1)
  sm(2) = c*sm(2)
  IF (x < 1.d-2) GO TO 200
  !                       TERMINATION
130 e = EXP(-sz(1))
  qf(1) = e*COS(-sz(2))
  qf(2) = e*SIN(-sz(2))
  ts(1) = qf(1)*sm(1) - qf(2)*sm(2)
  ts(2) = qf(1)*sm(2) + qf(2)*sm(1)
  sm(1) = ts(1)
  sm(2) = ts(2)
  IF (mo /= 0) GO TO 140
  w(1) = sn*(0.5d0 + (0.5d0 - sm(1)))
  w(2) = - sn*sm(2)
  RETURN
140 IF (sn == 1.d0) GO TO 141
  w(1) = 2.d0 - sm(1)
  w(2) = -sm(2)
  RETURN
141 w(1) = sm(1)
  w(2) = sm(2)
  RETURN
  !               MODIFIED ASYMPTOTIC EXPANSION
200 e = EXP(-sz(1))
  qf(1) = e*COS(-sz(2))
  qf(2) = e*SIN(-sz(2))
  w(1) = qf(1)*sm(1) - qf(2)*sm(2)
  w(2) = qf(1)*sm(2) + qf(2)*sm(1)
  IF (mo == 0) GO TO 210
  w(1) = 1.d0 + sn*w(1)
  w(2) = sn*w(2)
  RETURN
210 IF (sn < 0.0) RETURN
  w(1) = - w(1)
  w(2) = - w(2)
  RETURN
END SUBROUTINE dcerf



PURE FUNCTION dnorm(x, y) RESULT(fn_val)
  ! Replaces the statement function anorm in the F77 code.
  REAL (8), INTENT(IN) :: x, y
  REAL (8)             :: fn_val
  fn_val = MAX(ABS(x), ABS(y))
  RETURN
END FUNCTION dnorm



PURE SUBROUTINE dcrec (x, y, u, v)
  !-----------------------------------------------------------------------
  !             COMPLEX RECIPROCAL U + I*V = 1/(X + I*Y)
  !-----------------------------------------------------------------------
  REAL (8), INTENT(IN)   :: x, y
  REAL (8), INTENT(OUT)  :: u, v
  REAL (8) :: d, t
  IF (ABS(x) > ABS(y)) GO TO 10
  t = x/y
  d = y + t*x
  u = t/d
  v = -1.d0/d
  RETURN
10 t = y/x
  d = x + t*y
  u = 1.d0/d
  v = -t/d
  RETURN
END SUBROUTINE dcrec



PURE SUBROUTINE cdivid (ar, ai, br, bi, cr, ci)
  !-----------------------------------------------------------------------
  !     REAL (8) COMPLEX DIVISION C = A/B AVOIDING OVERFLOW
  !-----------------------------------------------------------------------
  REAL (8), INTENT(IN)   :: ar, ai, br, bi
  REAL (8), INTENT(OUT)  :: cr, ci
  REAL (8) :: d, t, u, v
  IF (ABS(br) <= ABS(bi)) GO TO 10
  t = bi/br
  d = br + t*bi
  u = (ar + ai*t)/d
  v = (ai - ar*t)/d
  cr = u
  ci = v
  RETURN
10 IF (bi == 0.d0) GO TO 20
  t = br/bi
  d = bi + t*br
  u = (ar*t + ai)/d
  v = (ai*t - ar)/d
  cr = u
  ci = v
  RETURN
  !     DIVISION BY ZERO. C = INFINITY
20 cr = HUGE(1.0d0)
  ci = cr
  RETURN
END SUBROUTINE cdivid






PURE SUBROUTINE erfcm2 (mo, z, w)
  !-----------------------------------------------------------------------
  !           CALCULATION OF ERFC(Z) USING THE TAYLOR SERIES
  !                          AROUND Z0 = 2
  !-----------------------------------------------------------------------
  INTEGER, INTENT(IN)     :: mo
  REAL (8), INTENT(IN)   :: z(2)
  REAL (8), INTENT(OUT)  :: w(2)
  REAL (8) :: eps, h(2), t(2), tol, x, y
  INTEGER   :: j, n
  !------------------------------
  !     C = (2/SQRT(PI))*EXP(-4)
  !     E = ERFC(2)
  !------------------------------
  REAL (8), PARAMETER :: c = .020666985354092053857068941306585476d0,  &
       e = .0046777349810472658379307436327470714d0
  !------------------------------
  REAL (8) :: a(63) 
  ! = (/ .20000000000000000000000000000000000D+01,  &
  !      .23333333333333333333333333333333333D+01,  .16666666666666666666666666666666667D+01,  &
  !      .63333333333333333333333333333333333D+00, -.22222222222222222222222222222222222D-01,  &
  !      -.16349206349206349206349206349206349D+00, -.76984126984126984126984126984126984D-01,  &
  !      -.24250440917107583774250440917107584D-02,  .12716049382716049382716049382716049D-01,  &
  !      .50208433541766875100208433541766875D-02, -.25305969750414194858639303083747528D-03,  &
  !      -.78593217482106370995259884148773038D-03, -.19118154038788959423880058800693721D-03,  &
  !      .46324144207742091339974937858535742D-04,  .33885549097189308829520469732109944D-04,  &
  !      .28637897646612243562134629672756034D-05, -.29071891082127275370004560446169188D-05,  &
  !      -.89674405786490646425523560263096103D-06,  .96069103941908684338469767911200105D-07,  &
  !      .99432863129093191401848891268744113D-07,  .97610310501460621303387795457283579D-08,  &
  !      -.65557500375673133822289344530892436D-08, -.18706782059105426900361744016236561D-08,  &
  !      .20329898993447386223176373714372370D-09,  .16941915827254374668448114614201210D-09,  &
  !      .10619149520827430973786114446699534D-10, -.10136148256511788733365237088810952D-10,  &
  !      -.21042890133669970575386166675721692D-11,  .37186985840699828780916522245407087D-12,  &
  !      .17921843632701679986488128324051002D-12, -.89823991804248069863542565948598397D-16,  &
  !      -.10533182313660970970232171410372199D-13, -.12340742690978398320850088252659714D-14,  &
  !      .44315624546581333350568244777175883D-15,  .11584041639989442481950487524296214D-15,  &
  !      -.10765703619385988116658460442219647D-16, -.70653158723054941879586082239984222D-17,  &
  !      -.18708903154917138727191793341667090D-18,  .32549879318817103966053527398133297D-18,  &
  !      .40654116689599228385911733319215613D-19, -.11250074516817311101947327325293424D-19,  &
  !      -.28923865378584966737386008432031980D-20,  .23653053641701517160704870522922706D-21,  &
  !      .14665384680061888088099002254334292D-21,  .26971039707314316218154193225264469D-23,  &
  !      -.58753834789274356433279671015522650D-23, -.59960357240498652932299485494869633D-24,  &
  !      .18586826578121663981412155416486531D-24,  .38364131854692721721867481914852428D-25,  &
  !      -.41342210492630142578080062451711039D-26, -.17646283105274988992381528904600860D-26,  &
  !      .19828685934364181151988692232131608D-28,  .65592252170840353572672782446212733D-28,  &
  !      .40626551379996340638338449938639730D-29, -.20097984104191034713653294173834095D-29,  &
  !      -.28104226475997460044096389060743131D-30,  .48705319298749358709127987806547949D-31,  &
  !      .12664655832830787747161769929972617D-31, -.75168312488894341862391776330113688D-33,  &
  !      -.45760473722605993842481669806804415D-33, -.56725491529575395930156379514718000D-35,  &
  !      .13932664042920082608489441616061541D-34,  .10452448992516358449586503951463322D-35 /)

  a = (/ .20000000000000000000000000000000000D+01,  &
       .23333333333333333333333333333333333D+01,  .16666666666666666666666666666666667D+01,  &
       .63333333333333333333333333333333333D+00, -.22222222222222222222222222222222222D-01,  &
       -.16349206349206349206349206349206349D+00, -.76984126984126984126984126984126984D-01,  &
       -.24250440917107583774250440917107584D-02,  .12716049382716049382716049382716049D-01,  &
       .50208433541766875100208433541766875D-02, -.25305969750414194858639303083747528D-03,  &
       -.78593217482106370995259884148773038D-03, -.19118154038788959423880058800693721D-03,  &
       .46324144207742091339974937858535742D-04,  .33885549097189308829520469732109944D-04,  &
       .28637897646612243562134629672756034D-05, -.29071891082127275370004560446169188D-05,  &
       -.89674405786490646425523560263096103D-06,  .96069103941908684338469767911200105D-07,  &
       .99432863129093191401848891268744113D-07,  .97610310501460621303387795457283579D-08,  &
       -.65557500375673133822289344530892436D-08, -.18706782059105426900361744016236561D-08,  &
       .20329898993447386223176373714372370D-09,  .16941915827254374668448114614201210D-09,  &
       .10619149520827430973786114446699534D-10, -.10136148256511788733365237088810952D-10,  &
       -.21042890133669970575386166675721692D-11,  .37186985840699828780916522245407087D-12,  &
       .17921843632701679986488128324051002D-12, -.89823991804248069863542565948598397D-16,  &
       -.10533182313660970970232171410372199D-13, -.12340742690978398320850088252659714D-14,  &
       .44315624546581333350568244777175883D-15,  .11584041639989442481950487524296214D-15,  &
       -.10765703619385988116658460442219647D-16, -.70653158723054941879586082239984222D-17,  &
       -.18708903154917138727191793341667090D-18,  .32549879318817103966053527398133297D-18,  &
       .40654116689599228385911733319215613D-19, -.11250074516817311101947327325293424D-19,  &
       -.28923865378584966737386008432031980D-20,  .23653053641701517160704870522922706D-21,  &
       .14665384680061888088099002254334292D-21,  .26971039707314316218154193225264469D-23,  &
       -.58753834789274356433279671015522650D-23, -.59960357240498652932299485494869633D-24,  &
       .18586826578121663981412155416486531D-24,  .38364131854692721721867481914852428D-25,  &
       -.41342210492630142578080062451711039D-26, -.17646283105274988992381528904600860D-26,  &
       .19828685934364181151988692232131608D-28,  .65592252170840353572672782446212733D-28,  &
       .40626551379996340638338449938639730D-29, -.20097984104191034713653294173834095D-29,  &
       -.28104226475997460044096389060743131D-30,  .48705319298749358709127987806547949D-31,  &
       .12664655832830787747161769929972617D-31, -.75168312488894341862391776330113688D-33,  &
       -.45760473722605993842481669806804415D-33, -.56725491529575395930156379514718000D-35,  &
       .13932664042920082608489441616061541D-34,  .10452448992516358449586503951463322D-35 /)

  !------------------------------
  !     ****** EPS IS A MACHINE DEPENDENT CONSTANT. EPS IS THE
  !            SMALLEST NUMBER SUCH THAT 1.d0 + EPS .GT. 1.d0
  eps = EPSILON(1.0d0)
  !------------------------------
  tol = eps*1.d+12
  h(1) = 1.d0 + (1.d0 - z(1))
  h(2) = - z(2)
  x = 1.d0
  y = 0.d0
  w(1) = a(30)
  w(2) = 0.d0
  DO n = 31,63
     t(1) = x*h(1) - y*h(2)
     t(2) = x*h(2) + y*h(1)
     x = t(1)
     y = t(2)
     t(1) = a(n)*x
     t(2) = a(n)*y
     w(1) = w(1) + t(1)
     w(2) = w(2) + t(2)
     IF (dnorm(t(1), t(2)) <= tol*dnorm(w(1), w(2))) EXIT
  END DO
  DO j = 1,29
     n = 30 - j
     x = h(1)*w(1) - h(2)*w(2)
     w(2) = h(1)*w(2) + h(2)*w(1)
     w(1) = a(n) + x
  END DO
  x = h(1)*w(1) - h(2)*w(2)
  w(2) = h(1)*w(2) + h(2)*w(1)
  w(1) = 1.d0 + x
  x = c*(h(1)*w(1) - h(2)*w(2))
  w(2) = c*(h(1)*w(2) + h(2)*w(1))
  w(1) = e + x
  IF (mo == 0) RETURN
  !                     COMPUTE EXP(Z*Z)*ERFC(Z)
  x = z(1)*z(1) - z(2)*z(2)
  y = 2.d0*z(1)*z(2)
  x = EXP(x)
  t(1) = x*COS(y)
  t(2) = x*SIN(y)
  x = t(1)*w(1) - t(2)*w(2)
  y = t(1)*w(2) + t(2)*w(1)
  w(1) = x
  w(2) = y
  RETURN
END SUBROUTINE erfcm2
